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The behaviors of the animals or embodied agents are characterized by the dynamic 
coupling between the brain, the body, and the environment. This implies that control, 
which is conventionally thought to be handled by the brain or a controller, can partially be 
outsourced to the physical body and the interaction with the environment. This idea has 
been demonstrated in a number of recently constructed robots, in particular from the field 
of "soft robotics" Soft robots are made of a soft material introducing high-dimensionality, 
non-linearity, and elasticity, which often makes the robots difficult to control. Biological 
systems such as the octopus are mastering their complex bodies in highly sophisticated 
manners by capitalizing on their body dynamics. We will demonstrate that the structure 
of the octopus arm cannot only be exploited for generating behavior but also, in a sense, 
as a computational resource. By using a soft robotic arm inspired by the octopus we 
show in a number of experiments how control is partially incorporated into the physical 
arm's dynamics and how the arm's dynamics can be exploited to approximate non-linear 
dynamical systems and embed non-linear limit cycles. Future application scenarios as well 
as the implications of the results for the octopus biology are also discussed. 
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1. INTRODUCTION 

Biological systems have certain morphologies' and material char- 
acteristics that improve their adaptivity and increase their proba- 
bility to survive. This suggests that control is not only located in 
the brain, but that there is a tight coupling between the brain, 
the body, and the environment, an idea that is usually termed 
embodiment (Pfeifer and Bongard, 2006; Pfeifer et al, 2007, 2012; 
Li et al, 2011b; Nakajima et al, 2011c, 2012a,b). Recently moti- 
vated by the fact that soft material is ubiquitous in the body 
structures of living creatures, a new family of robots, soft robots, 
has been constructed with the aim of incorporating flexible ele- 
ments (Trivedi et al., 2008; Steltz et al., 2009; Brown et al, 2010; 
Shepherd et al., 2011; Pfeifer et al., 2012). These robots have sig- 
nificant advantages over traditional articulated robots in terms of 
morphological flexibility and interactional safety (Trivedi et al., 
2008; Li et al., 2011b). However, controlling them with conven- 
tional techniques is difficult because of their high-dimensional 
body structures and their diverse body dynamics, which are due to 
their non-linearity and elasticity. In this context, the octopus has 
been a good source of inspiration for roboticists to learn a control 
strategy for soft robots. An octopus has hyper-redundant limbs 
with a virtually unlimited number of degrees of freedom (DOT), 
and its movements are known to be highly sophisticated (Sumbre 



By morphology, we do not only refer to the shape, but also sensor and 
actuator distributions, and physical properties, such as stiffness, etc. 



et al., 2001, 2005; Trivedi et al, 2008). From a conventional 
control perspective, the octopus's method of controlling move- 
ment and harnessing its non-linear body dynamics is outstanding 
and far-reaching. 

It is well known that the nervous system of the octopus is 
highly distributed throughout the entire body. It has a relatively 
small central brain (about 50 million neurons), a central nervous 
system (CNS), which controls the large peripheral nervous sys- 
tem (PNS) of the arms (about 300 million neurons), integrates 
information from the visual system, and then issues commands 
to lower motor centers controlling the elaborate neuromuscular 
system of the arms. A typical example showing the effectiveness 
of this distribution of the nervous system is the reaching behav- 
ior (Gutfreund et al, 1996; Gutfreund, 1998; Sumbre et al, 2001; 
Yekutieli et al, 2005a,b). Reaching behavior consists of a bend 
propagation along the arm toward the tip in a highly stereotypical 
and invariant way. Sumbre et al. showed that the arm extensions 
can be evoked in arms whose connection with the brain have 
been severed (Sumbre et al., 2001). Because the evoked motions 
in denervated octopus arms were qualitatively and kinematicaUy 
identical to natural bend propagations, an underlying motor pro- 
gram appears to be embedded in the neuromuscular system of 
the arm, which does not require continuous central control (Li 
et al., 201 la, 2012, 2013; Nakajima et al., 2011a,b; Kuwabara et al, 
2012). In addition, the muscular organization of the octopus's 
arm has a characteristic structure called muscular-hydrostats (Kier 
and Smith, 1985; Smith and Kier, 1989; Taylor and Kier, 2003; 
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Feinstein et al., 201 1). In such structures, the volume of the organ 
remains constant during all movements, enabling the muscles 
themselves to perform all the functions usually performed by the 
skeleton (Sumbre et al, 2001, 2005; Taylor and Kier, 2003). This 
suggests that the body of the octopus arm is highly involved in the 
production of movements. Accordingly, in robotics, there have 
been several attempts to characterize the role of the muscular- 
hydrostat system in terms of its anatomical structure (Mazzolai 
et al, 2007; Laschi et al, 2009, 2012; Vavourakis et al, 2012a,b) 
and functionality (Nakajima et al., 2013). 

In this paper, along the lines of these biological findings, we 
aim to provide one quantitative evidence that the structure of the 
octopus's arm has the potential to embed multiple motor pro- 
grams without any support from the external controller. Recently, 
it has been shown that non-linear mass-spring networks can be 
used as a computational resource (Caluwaerts and Schrauwen, 
2011; Hauser et al, 2011, 2012; Sumioka et al., 2011; Caluwaerts 
et al., 2013; Nakajima et al., 2013). These works imply that the 
non-linear and elastic body dynamics of soft robots are not 
drawbacks for control, but rather can be directly exploited as a 
computational resource. In this paper, we build on theoretical 
models (Hauser et al., 2011, 2012) that have been proposed in 
the context of reservoir computing. 

The term reservoir computing has been proposed by 
Schrauwen et al. (2007) for a set of machine learning tech- 
niques used to emulate complex, non-linear computations. The 
idea is to drive a high-dimensional, non-linear dynamical system 
(which has been randomly initialized, but afterwards fixed) with 
a low-dimensional input stream. This dynamical system, typically 
referred to as the reservoir, provides highly complex, but repro- 
ducible responses in its state space to the input. It operates as a 
type of temporal and finite "kernel" by projecting non-linearily 
the low-dimensional input into the high-dimensional state space 
of the reservoir. Furthermore, since a reservoir consists of dynam- 
ical systems, it exhibits a memory, which fades out exponentially 
(i.e., fading memory)^. A remarkable property of the approach 
is, if the reservoir is complex enough (i.e., high-dimensional 
and non-linear), it can been shown that it is sufficient to add 
a simple linear, static readout from the high-dimensional state 
space to emulate non-linear complex computations. Such reser- 
voir computing setups have been proven to outperform other 
machine learning techniques in a number of difficult tasks; see 
Jaeger (2003) for example. Another remarkable property of this 
setup is that the requested properties for computationally pow- 
erful reservoirs turn out to be rather general. Hence, a number 
of different implementations for reservoirs have been proposed 
(Schrauwen et al., 2007). For example, simple, abstract dynamical 
systems are used for echo state networks (Jaeger, 2002; Verstraeten 
et al., 2007; Lukosevicius and Jaeger, 2009), or models of neurons 
are used in liquid state machines (Maass et al., 2002). Lately, it 
has been demonstrated that complex, compliant bodies of bio- 
logical systems and robots have the potential to serve as such 
a reservoir as well, see Hauser et al. (2011) and Hauser et al. 
(2012). 



The memory is due to the integration capability of dynamical systems (i.e., 
It accumulates information over time). 



Here, we demonstrate that the soft robotic arm inspired by the 
octopus can be used as a reservoir. This means, by simply attach- 
ing a static, linear readout from the high-dimensional non-linear 
dynamics of the octopus arm, one can emulate complex, non- 
linear computation without altering the physical system itself 
That is, we employ the physical body as part of a computational 
device. In this paper, a 3D dynamic model of this soft robotic 
arm is used as a test platform. Compared with the model used 
in Hauser et al. (2011) and Hauser et al. (2012), our model is 
more biologically inspired and physically feasible. It is a mass- 
spring-damper system, where the springs are aligned to emulate 
the octopus's muscular organization, and embeds the characteris- 
tic properties of a muscular-hydrostat. The arm is also assumed 
to be immersed in an underwater environment, in which the 
water friction constants are approximated by the computational 
fluid dynamics (CFD) simulations. As a result, the arm reveals 
highly non-linear body dynamics when actuated. By using this 
platform, we demonstrate that its body dynamics can be exploited 
as a computational resource. To test its power, we defined two 
types of tasks: first, to emulate complex non-linear dynamical 
systems, where we investigate whether the body dynamics are 
exploited as a computational resource; second, to implement a 
closed-loop control. We used several non-linear limit cycles to 
see how they can be embedded directly into the soft robotic arm 
without any support from an external controller. The choice of 
example functions adopted for each type of task is motivated to 
evaluate the non-linearity and memory that the body dynamics 
contains. 

This paper is organized as follows. In section 2, we start 
by explaining the overall setting of the 3D dynamic model of 
the soft robotic arm platform and show how the arm emulates 
the muscular organization of the octopus in detaU. The input- 
output relations adopted and the experimental procedures are 
also provided in detail. In section 3, we explain the results for 
a series of tasks in detail, and in section 4, we give concluding 
remarks, including future extension scenarios of our proposed 
approach. 

2. MATERIALS AND METHODS 

In this section, we provide a detailed description of the soft 
robotic arm simulator model and explain how to exploit the sys- 
tem as a computational resource by introducing input-output 
relations. The experimental procedure is also explained in detail. 

2.1. DYNAMIC MODEL OF A SOFT ROBOTIC ARM INSPIRED BY THE 
OCTOPUS 

In this paper, we use a 3D dynamic model of a soft robotic arm 
inspired by the octopus (Kang et al, 2011, 2012). The model is 
currently applied for testing purposes for control architectures of 
soft robotic arms (Kuwabara et al., 2012; Nakajima et al, 2012a), 
and has been validated to have good agreement with a physical 
soft robotic arm platform (Zheng et al, 2012). The overall struc- 
ture of the entire arm is shown in Figure IF. It is assumed to be 
immersed in an underwater environment, and the base of the arm 
is able to rotate in any direction. It consists of 20 compartments, 
and each compartment implements the unique characteristics of 
octopus muscles, called muscular-hydrostats. In an octopus arm. 
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FIGURE 1 I (A) Histological transverse section of an octopus arm, showing 
the longitudinal muscles (L), the transverse muscles (T), the nerve fibers (N), 
and the oblique muscles (O). (B) Overall structure of the muscular-hydrostat 
system (a single compartment) used in this paper. (C) Schematics showing a 
longitudinal spring, a transverse spring (D), and a ceiling plane (E). In (C), 
fdamp = Ciifli, 4,iff = Kiijili - l^..). and i;,- = 5;^,- + ^vi,' In (D), fdamp = Cru'l^-, 



4,ijf = /(njC/J — /g,y), and S|y = 8„y + 8„y. See the text for details. (F) The entire 
soft robotic arm inspired by the octopus. It consists of 20 compartments, 
which are numbered from 1 to 20 from the base to the tip. Each 
compartment contains four longitudinal springs and one transverse disk. The 
blue line represents the line connecting the centers of each transverse disk. 
The base of the arm is able to rotate in any direction. See the text for details. 



muscles are organized into transverse, longitudinal, and obliquely 
oriented groups (Figure lA). This special muscular organization 
forms the structures of the muscular-hydrostats. Their main 
property is that their volume remains constant during muscle 
contractions. The result is that if the diameter of the muscular- 
hydrostats decreases, then their length increases, and vice versa. 
Several proposed models deal with the muscular-hydrostat sys- 
tem of the octopus [e.g.. See Yekutieli et al. (2005a,b) and Kang 
et al. (2012)]. The overall structure of the muscular-hydrostat 



system adopted in this paper is shown in Figure IB. We begin our 
description by focusing on the model of a single compartment, 
and then progress to describing an entire arm. 

A single compartment is a mass-spring-damper system, 
shaped like a circular truncated cone, consisting of a base plane, 
a ceiling plane with four transverse springs, a central strut, and 
four longitudinal springs, which emulate the anatomical struc- 
ture of the muscle alignment in a real octopus arm. (Note that, 
although we use the term "spring," it is a model for a muscle. 
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SO it has mass and damping.) The longitudinal springs con- 
trol the position and orientation of the ceiling plane, while 
the transverse springs control the radius of the ceiling plane. 
The central strut provides kinematic constraints to guarantee 
the unique position of the ceiling plane. It is considered as 
an ideal prismatic joint without mass, damping, and stiffness. 
The system has an isovolmetric structure, which provides forces 
constantly aiming to maintain its volume and is an expression 
of the property of the muscular-hydrostats, and thus, all the 
springs are assumed to be implicitly or explicitly coupled with 
each other. The values for all the parameters of the model (e.g., 
spring coefficients, damping, etc.) are either inspired by the octo- 
pus or directly drawn from biological data (Kier and Curtin, 
2002; Lieber, 2002; Shinohara et al, 2010). Standard SI units 
are used for the variables, and all of the ordinal differential 
equations presented are solved using the 4th order Runge-Kutta 
method, where dt is set to 0.001 s for the system throughout this 
paper. 

Coordinates are defined on a base plane and a ceiling 
plane, denoted by Aj(xj, yj, zj) and Bj{uj, vj, wj), respectively 
(Figure IB), where is the index number of the compartment. 
A vector expressing a longitudinal spring (d;,) is: dy = pj -|- by — 
ay, where py = AjBj = [0 0 hj]^ is the position vector of the cen- 
ter of the ceiling plane, by is the position vector of joint By, 
and ay = AjAy is the position vector of joint Ay. The length of 
the ith spring in compartment j is I'-j and can then be obtained 

by l'^ = ^djj ■ dy. The dynamics of a longitudinal spring is 
expressed by: 



spring. The dynamics of the length of the transverse spring are 
described as: 



+ K, 



(4) 



where h^j is the control force, 8,,y is the isovolumetric constraint 
force (which will be discussed in detail later), Fcli) = — FfiLy is 
the joint force of By acting on the ceiling plane, Uf,y is the unit 
vector of by, M^y is the mass of the spring, V^^^ is the initial radius 
of the ceiling plane, C,-,-,- is the damping coefficient, and K„j is the 
stiffness coefficient (Figure ID). Figure IE shows the illustration 
of the ceiling plane. The equation for the motion of the ceiling 
plane is: 



Mce,0 [O 0 h^^ = FcLl; + FcL2; + FcLS; 

+ FcL4j + FCQ + 



(5) 



where hj is the height of the ceiling plane of compartment 
FcQ = — Fbc; is the joint force on Bj acting on the ceUing plane, 
and Fexj is the external force. The rotation is formulated as: 



''^Iceil; [Pi 0]^ = ^ J5,by X ^iRj.FcL,j 



ex; 5 



(6) 



where %cij is the control force, is the isovolumetric constraint 
force, which will be also explained further in Equation (7),/BLzy is 
the component force of joint 5y acting along the spring i of com- 
partment;', is the initial length of the spring, M;y is the mass of 
the spring, C/y is the damping coefficient, and Knj is the stiffness 
coefficient (Figure IC). Then, the rotation of the longitudinal 
springs can be formulated in frame A by: 



where ^'Tccj is the constraint torque of joint Bj acting on the ceil- 
+ Kiij — Zq-^ , (1) ing plane, is the Euler rotation matrix, ^^by is the position 

vector of By expressed in frame Bj, and ^^ Iceilj is the inertia matrix 
of the ceiling plane. 

As explained previously, the system is isovolumetric due to 
its muscular hydrostat structure. This means that an increase in 
the longitudinal length will result in a reduction in the cross- 
sectional area and vice versa. A pair of antagonistic forces are 
applied to the longitudinal and transverse springs to guarantee 
the isovolumetric constraints. These are expressed as: 



7yM;y 



(2) 



where f^ixyij is the component force of joint By acting perpendic- 
ular to the spring, 7y is the inertia moment of the spring about 
Ay, and co/y is the angular velocity of the spring about Ay , where 
mij = (dij/llj) X (vBy//y), and VBy is the velocity of By. 

To interlink several compartments serially the reaction forces 
between the longitudinal springs and the base, F^^y, need to be 
calculated. These reaction forces are obtained by: 



FaIi; = Mujdij - Fbl,;, 



(3) 



where ^AUj is the joint force on Ay, and Fbl,; = i^Lxyij + fBLzy is 
the joint force on By . 

The four transverse springs spread around the central point of 
the ceiling plane. Figure ID shows the illustration of a transverse 



-Ki,x |8,y| X (V,j-Vo,), 



-Krv X 



+ F„j • Pj/hj 



(7) 



x(Vg-\^Oj), (8) 



where Vq is the actual volume of compartment j, Vgj is the 
initial volume of the compartment, Kiy is the constraint force 
gain for the longitudinal springs, and 7C„ is the constraint force 
gain for the transverse springs. From Equation (7), it can be seen 
that the constraint force is a function of the transverse spring 
force, hcij, and the compartment volume change, Vq — Vq;- By 
applying %„j to Equation (1), the longitudinal springs will act 
against the transverse springs to drive the volume change to zero. 
Similarly, another constraint force 8,,/, is obtained by Equation 
(8) and applied to Equation (4) for the transverse springs to can- 
cel the volume change induced by the longitudinal springs. Note 
that the external force, Fexj, is included in Equation (8) because 
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it is acting on the compartment on joint Bj, which may also 
change the length of the longitudinal springs and the volume of 
the compartment. 

In addition to the forces generated by the muscles, typical 
external forces applied to the soft robotic arm in an underwa- 
ter environment are gravity, buoyancy, and hydrodynamic forces. 
These are considered as distributed forces acting on each com- 
partment as: 

= ¥gj + ¥bj + (FhydDj + FhydLj) . (9) 

where Fexj is the total external force acting on compartment j, ¥gj 
is the gravity force, F^j is the buoyancy force, and V^jdj is the 
hydrodynamic force composed of the water drag force, FhydDi> 
and the water lift force, FhydLj- 

The direction of buoyancy always opposes gravity. Thus, the 
resulting force due to gravity and buoyancy is: 

¥gj + = (Po - PwWcjgUgj = PoVggeUgj, (10) 

where p„, is the density of water, pg is the density of the octo- 
pus arm, Vg is the volume of compartment j, Ugj is the unit 
vector indicating the direction of gravity for compartment ;', g 
is the gravity constant, and gi- is the equivalent gravity constant 
defined as: 

By adjusting the value of gf,, both gravity and buoyancy forces are 
included in the model. 

The hydrodynamic forces applied to an octopus arm during 
a movement through a fluid medium are shown in Figure 2. For 
compartment j, the drag force FhydDj is parallel to the velocity vec- 
tor Vj of the fluid (or the uniform arm velocity in a stationary 
fluid) and the lift force FhydL; is perpendicular to V,-, according to 

1 2 
FhydDj = -CojArjPwWVjW Uyj, (12) 

FhydLj = ^CljArjpJlVjW^ujj, (13) 

where Cdj and Cij are the drag and lift coefficients, respectively, 
Arj is the reference area of compartment j, Uyj is the unit vector 
indicating the direction of V,, and is the unit vector of nor- 
mal direction for Vj. The hydrodynamic force coefficients, Coj 
and Cij, for a segmented arm are obtained from high fidelity 
CFD simulations (Kazakidi et al., 2012). They were found to be 
dependent on the flow incidence angle Qj, and the configuration 
of the arm (e.g., straight vs. bended). As a first approximation of 
the hydrodynamic forces common in robotic literature (Ijspeert, 
2001; Kazakidi et al., 2012), dependence on arm configuration 
was ignored. Therefore, a single value for each coefficient at 
specific angles of Qj for a straight arm was identified by CFD 
simulations and approximated by a 4th order polynomial in the 



simulator, expressed as follows: 

Cdj = ePQl + ePe] + ePej + e^Qj + ef, ( 14) 

Clj = e[e^ + + 4eJ e^e^ 4, (15) 

where ej^ ^ and e^ ^ are the parameters identified by CFD simula- 
tions. The hydrodynamic forces for each compartment were then 
computed according to Equations (12) and (13), where Vj was 
taken as the velocity of Bj. All the parameters used in this study 
are shown in Table 1. 

This model is intrinsically non-linear. The non-linearities 
of the system are partly introduced by its kinematics (Kang 




FIGURE 2 I Hydrodynamic forces acting on the soft robotic arm. 



Table 1 | Parameters for the soft robotic arm adopted in this paper. 



Parameter 


Value 


Parameter 


Value 


'o,y (n^) 


0.0200 


/C„y (N/m) 


196 


loi, (n^) 


0.015-0.0067 


Krij (N/m) 


1570 


Miij (kg) 


1.25 x10-3 


(l/m3) 


1.0 xlO" 


Mrij (kg) 


1.25 x10-3 


Krv (l/m^) 


1.0 xlO" 


Qij (N s/m) 


1.0 


Voj (m3) 


1.38 x10-6 








-2.99 X 10-6 


C„y (N s/m) 


1.0 


Pw (kg/m^) 


1000 




-5.5 X 10-9 


^^ 


1.8 x10-3 


pD 


6.3 xlO-^ 


pi- 


-5.9 X 10-^ 


pD 


-7.7 X 10-6 




2.8 x10-5 


"a 


0.0015 


^4 


0.0011 


pD 


0.017 




0.00089 



Note that, for l^-j and Voj, they decrease monotonically, according to the 
compartment number in the given range in the table. 
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et al., 2012). The relation between the spring length and the 
ceiling plane posture (position and orientation) is non-linear. 
Therefore, the system dynamics become non-linear as well. Also, 
the calculation of the isovolumetric forces and hydrodynamic 
forces introduces non-linearities. See Kang et al. (2012) for a 
detailed discussion of the model. In the majority of cases, these 
non-linearities are undesirable from the viewpoint of classical 
control theory. However, as previously mentioned in section 1, 
such a complex body could potentially be used as part of a 
computational device, if appropriate inputs and readouts are 
applied. 

2.2. EXPERIMENTAL PROCEDURE 

Our aim in this paper is to demonstrate whether a soft robotic 
arm can be exploited as a computational resource, as well as a 
controller. Accordingly, we need to define inputs (In(f)) to the 
system and how to generate corresponding outputs (0(f -|- 1)). 
In this paper, we apply the position control of the base rota- 
tion as an input, and an output is generated by the weighted 
sum of the longitudinal spring lengths of all 20 compartments 
(Figures). 

Based on this I/O scheme, we set two types of tasks for 
our demonstrations. First, we consider the emulation tasks of 
non-linear dynamical systems (Figure 3A), which aims to show 
whether the soft robotic arm can be exploited as a computa- 
tional resource, including sufficient non-linearity and memory. 
The second task is to embed closed-loop control onto the soft 
robotic arm itself (Figure 3B). In particular, we aim to embed 
non-linear limit cycles, which are especially appealing for the 
control of robots. Typically, such limit cycles are implemented 



using non-linear oscillators, such as central pattern genera- 
tors (CPGs), or a network of such oscillators (Righetti and 
Ijspeert, 2008). We here aim to demonstrate that the body of 
the soft robotic arm itself can be used to generate such limit 
cycles. 

As explained in section 1, our approach is comparable to a 
reservoir computing approach, which normally uses randomly 
coupled non-linear elements as a computational resource (Jaeger, 
2002; Maass et al, 2002; Jaeger and Haas, 2004). In the conven- 
tional reservoir computing approach, since each computational 
element is coupled randomly, each element possess a uniform 
role in the computation in the statistical sense. On the other 
hand, if we exploit the robot's body as a reservoir, accord- 
ing to the intrinsic structure of the body, each part of the 
body shows qualitatively different dynamics, which may lead 
to specific role distributions corresponding to each body part. 
Accordingly, in this paper, we will investigate how the compu- 
tational role is distributed through the arm in each task. In 
the following subsections, we provide detailed descriptions for 
each task. 

2.2. 1. Task 1: non-linear dynamical system emulation tasks 

In order to evaluate the computational power of the system, we 
here set non-linear dynamical system emulation tasks, which are 
often used as benchmark tasks (Jaeger, 2002; Verstraeten et al, 
2007; Hauser et al, 2011) in the context of recurrent neural 
network learning (Atiya and Parlos, 2000) and the reservoir com- 
puting approach (Jaeger, 2002; Maass et al, 2002; Jaeger and 
Haas, 2004). Each task requires a certain degree of non-linearity 
and memory to be performed by the system. As explained above. 



A Emulation Tasks 

Input Output 
(ln(t)) ..----^ (0(t+1)) 

O <i 

base rotation 



i /■ linear readout 



cps\x 




B Closed-loop Control 

Outputi ^ , Output2 
1/^ /t i)\ base rotation ,^ IT .„ 

{Oi(t+1)) ,-■ -N {02(t+1)) 



1 !• ; 



linear efi^out 

■' ' <3>/ 



FIGURE 3 I Schematics explaining the tasks adopted in this paper. 

(A) Schematics showing the nonlinear dynamical system emulation tasks 
(Task 1). An input (In(t)) is projected as a base rotation angle, and 
accordingly, the soft robotic arm shows passive body dynamics. By 
setting a linear readout for each longitudinal sphng length in each 
compartment, an output (O(t-I-I)) is calculated as a weighted sum of all 
the spring lengths. By adjusting only the linear readout, we demonstrate 



whether the system can emulate complex nonlinear dynamical systems. 
(B) Schematics showing the implementation of closed-loop control 
(Task 2). Similar to Task 1, by adjusting the linear readouts, two variables 
of a nonlinear limit cycle are emulated and fed back as the next input to 
the system to generate the base rotation movement. Accordingly, the 
nonlinear limit cycle is embedded onto the arm in a closed-loop manner. 
See the text for details. 
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we first apply position control of the base rotation as an input 
(Figure 3A), expressed as follows: 



where y(t) denotes the output of the system. The third task is a 
discrete Volterra series, expressed as follows: 



e(f) = (t)(f) = Scale X In(f), (16) 
In(f) = 0.2 sin(2Tt/i f x dt) sin(2Tt/2f x dt) sin(2it/3f x dt), (17) 

where 9(f) and <\>{t) are base rotations at timestep t along the x- 
axis and y-axis, respectively. The parameter Scale linearly scales 
the raw input, In(f), to the specific range of the base angle 
[degree], —R < 9(f), (t>(f) < R- This scaling parameter can be 
freely chosen, but should be fixed throughout the experiment. 
The detailed setting will be explained in the section 3. The param- 
eters /i, /2, and fi are set to 2.11, 3.73, and 4.33, respectively. 
Similar inputs were adopted in Hauser et al. (20 11 ), Sumioka et al. 
(2011), and Nakajima et al. (2013). 

According to the base rotation, the arm generates passive body 
dynamics. The output of the system is calculated by using the 
resulting spring dynamics as follows: 



S,,(f), 



(18) 



where s/,(f) is the length of the longitudinal spring i (i = 1, 2, 3, 
and 4) in compartment j {j = 1,2, ... 20) at timestep f. When 
the input is projected as the base rotation at timestep f , the corre- 
sponding length of the spring, li, is collected to Sy (f). The linear 

readout weight, w'^^^, corresponds to each spring length. Overall, 
the dynamics of 80 (= 4 x 20) spring lengths are the expression 
of the body dynamics in this paper. 

In order to achieve the required computation, we only 
train the linear readout {Wg^^). Since we have 80 nodes, 

Su(t),S2i{t),S3i(t), Suit), Suit), , 5420(f), for the lengths of 

the spring at timestep f, by collecting the lengths of the 
springs for M timesteps, we can generate an 80 xM matrbc 
L. We also collect the corresponding target outputs for M 
in a matrix T. Then, the optimal readout weights, W = 

[Woif w^it'M'out'M'oJit'W^ut, •■•■'i^ou?]^ can be obtained by 
W = L*T, where L* is the Moore-Penrose pseudo-inverse, since 
L is not a square matrix in general. 

According to the sent input, the system should emulate the 
following three non-linear dynamical systems as outputs. We pre- 
pared three corresponding output nodes to the system whose 
linear readouts are trained separately for each task. (This pro- 
cedure is often called multitasking.) The first task is a 2nd order 
non-linear dynamical system, expressed as follows: 

y(t +l) = 0.4y(t) + 0.4y(f)y(f - 1) + 0.6ln^(t) + 0.1, (19) 

where y(t) denotes the output of the system. The second task is a 
10th order non-linear dynamical system, expressed as follows: 



i-F0.05y(f) 



y{t +1) = 0.3y(f) + 0.05y(f) [l^j'^'^ " I (20) 



r(f+i) 



200 200 

11=0X2 



(ti,t2)In(f-ti)In(f-T2), (22) 



0 



h(xi, T2) = exp 



(ti X dt — (Xi)^ (t2 X dt — 1x2)^ 



2aj 



+ 



2al 



(23) 



-M.5In(f-9)In(f)-|-0.1, 



(21) 



where A is a scaling parameter set to 0.0001, y(t) and h(xi, T2) 
denote the output of the system and a Gaussian kernel, respec- 
tively. The parameters, [L\, (jl2. o\, and 02 are set as M-i = (JL2 = 
0. 1 and 01 = 02 = 0.05. Any computational model that can emu- 
late the above dynamical systems should have a certain degree of 
memory and non-linearity. Simply put, emulation of a 10th order 
system requires more memory and non-linearity than a 2nd order 
system, and emulation of the Volterra task requires more than the 
10th order system. 

For the experimental procedure, the soft robotic arm is first 
set in the resting state with 9(f) = ^{t) = 0, and before begin- 
ning the experiment, we start to run the arm with Equation (16) 
for Tini timesteps. This phase is to set the different initial posi- 
tions of the arm for each experimental trial; Tini is randomly 
determined from 0 to 1000 timesteps for each trial. The actual 
experimental trial consists of 16,000 timesteps, where the first 
1000 timesteps are for washout, the following 10,000 timesteps 
are for the training phase, and the final 5000 timesteps are for 
the evaluation phase. After Tini timesteps, we continue run- 
ning the arm with Equation (16) and the actual experiment 
begins. By collecting the lengths of the spring and the corre- 
sponding target outputs for each task in the training phase, we 
train the linear readouts for three outputs by adopting the previ- 
ously explained procedure. By using the trained linear readouts, 
we evaluate the performance of the system output by calculat- 
ing the mean squared error (MSE), MSE = i Ef=i(0(f + 1) - 
y(t + 1))^, where n = 5000. We here compare the performance 
of the system with outputs generated by simple linear regres- 
sion, 0(f -|- 1) = a X In(f) -|- b, where a and b are trained by 
using the same time series as in the training phase. As is clear 
from the equation, since the linear regressor only uses the input 
to generate the output, which does not contain non-linearity 
and memory, any task performance of the system better than 
the linear regressor can be said that the required non-linearity 
and memory to perform the task is positively exploited from 
the system. 

2.2.2. Task 2: closed-loop control — embedding non-linear limit 
cycles 

As previously explained, in this task, we aim to embed non-linear 
limit cycles in a closed-loop manner. The major difference from 
Task 1 is that the outputs generated by the system are fed back 
to the system itself as a motor command (an input) for the next 
timestep (Figure 3B). In particular, as will be explained later, 
we here aim to embed several limit cycles, which each have two 
variables. Accordingly, the outputs generated for the next motor 
commands (namely, Ini(f) and In2(f) for each variable) are 
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projected to 6(f) and (t)(t), respectively. The situation is expressed 
as follows: 



9(0 = 
^(0 = 

Ini(f) 
In2(f) 



lni(0, 
In2(f), 

= Oi(f), 

= 02(0, 



Oi(f+l) 

02(f+l) 



E20 v^4 
j= 1 Z^i = 

E20 v^4 
;"= 1 2^1 = 



,(f), 



(24) 



(25) 



(26) 



where w'^^^ j and w'^^^ ^ the linear readouts corresponding 
to the two outputs, Oi(f) and O2(0> respectively. (Note that, in 
Equation (24), unlike Equation (16) in Task 1, the scaling param- 
eter Scale is not introduced. As will be explained later, this is 
because we here aim to emulate the limit cycles, which are already 
scaled with a certain parameter value. So the scaling procedure is 
already included in the target outputs.) 

As in the procedure explained in Task 1, to train the system, we 
only adjust the linear readouts, w^^j j and w^^^ 2- However, the 
procedures differ in two points; first, during the training phase, 
we clamp the feedbacks from the system outputs, and provide the 
target outputs as inputs, which means, in Equation (25), we set 
Ini(f) = x\{t), In2(f) = X2(f), where x\{t) and X2(f) are the tar- 
get outputs at timestep f . Thus, the training phase is carried out 
with an open-loop, where the system was forced into the desired 
operative state by the target signals (this is called, teacher jorcin^ 
(Hauser et al., 2012). Second, when collecting the lengths of the 
springs in the training phase, we add white noise in the range of 
[—V, v]. By doing this, we can expect that the obtained optimal 
readouts will generate the target outputs even under the influ- 
ence of noise (Hauser et al., 2012). The appropriate degree of v is 
determined heuristically for each task. 

Here, we aim to embed three non-linear limit cycles. The 
first one is the dynamical systems of the Van der Pol equations, 
expressed as follows: 



Xi = X2, 

xi = -xi -I- (l -xfjxj. 



(27) 



The second one is a limit cycle, which we call the quadratic limit 
cycle (Hauser et al., 2012; Khalil, 2002), expressed as follows: 



X\ = Xl + X2 — 5X\ (X[ -|- , 
X2 = -2Xl +X2-X2{xj+ X^) . 



(28) 



The third one is a Lissajous curve with a frequency ratio of 
/1//2 = 2, expressed as follows: 



Xl = sin(fif), 
X2 = sin(/2f). 



(29) 



Since each limit cycle is symmetric about the point (0, 0) , we select 
the variable with the larger range, and scale both variables to the 



desired range of the base rotation [degree], — -R < Q(t),^(t) < 
R. As in Task 1, this scaling parameter can be freely chosen, 
but should be fixed throughout the experiment. Thus, what 
the system should emulate here is Xj(f -|- 1) = Scale x xi(t + 1) 
and ^(f + 1) = Scale x X2(t + 1) as Oi(f-|- 1) and 02{t + 1), 
respectively. Further settings on the parameter Scale will be 
explained in the section 3. For the Van der Pol system and the 
quadratic limit cycle, the ordinal differential equations are solved 
for each simulation timestep by using the 4th order Runge-Kutta 
method, where dt is set to 0.01. Note that the timescale of the arm 
model and these limit cycle is different. When we refer to time s, 
we always fixed our expression to the timescale to the arm model, 
otherwise we use the expression of simulation timestep to avoid 
the confusion. 

In the experimental procedure, the soft robotic arm is first 
set in the resting state with 6(f) = (t)(f) = 0, as in Task 1. We 
run the system with the teacher forcing condition for 70,000 
timesteps, and by discarding the first 10,000 timesteps, we use 
60,000 timesteps as for the training phase with white noise of 
degree v added in the spring lengths. After obtaining the opti- 
mal readout weights from these collected data, we initialize the 
arm's state to the resting state and again start to run the system 
with the teacher forcing condition. After 5000 timesteps of run- 
ning, we switch the inputs to the system output generated by the 
trained readout weights (Equations 25 and 26) and check whether 
it could embed the target limit cycle. Unlike Task 1, multitasking 
cannot be adopted (due to the feedback control), so each limit 
cycle is trained separately as a different trial. 

3. RESULTS 

In this section, we present the results of each task applied to 
indicate the performance of our system. We would like to note 
again that all the tasks presented in this section is performed with 
"one body", the same soft robotic arm explained in section 2.1, 
where all the parameter settings of the arm is fixed throughout 
the experiments. In addition, for Task 1, emulations of three non- 
linear dynamical systems are simultaneously performed for each 
experimental run (i.e., multitasking). 

3.1. TASK 1: NON-LINEAR DYNAMICAL SYSTEM EMULATION TASKS 

Figure 4 shows a typical example of the time series of In(f), the 
lengths of the springs, and the performance of each task during 
the evaluation phase. The plots show the case when R is set to 60. 
We can clearly see that, according to the input projected to the 
base rotation, our soft robotic arm shows diverse passive body 
dynamics (Figures 4A,B). Regarding the task performance, the 
system output shows better performances than the linear regres- 
sor in all the tasks (Figure 4C). In the emulation task for the 2nd 
order system, the linear regressor also showed relatively good per- 
formance. However, as we can see in the plots, as the degree of 
non-linearity and memory of the task increases in the 10th order 
system and Volterra task, the performance of the linear regressor 
decreases significantly (Figure 4C). On the other hand, our sys- 
tem shows relatively good performance even in the Voletrra series 
emulation task (Figure 4C). Table 2 shows the statistical compar- 
isons between the MSE of the output of the system and that of 
the linear regressor for each task. The values show the averaged 
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FIGURE 4 I A typical example of the task performance for Task 1 in 
terms of time series in the evaluation phase. (A) The time series of In(f). 
(B) Corresponding body dynamics, which are expressed as the time series 
of the lengths of springs 1, 2, 3, and 4. The spring dynamics of all 20 
compartments are overlaid. (C) Comparisons of the performance of the 
system output and the linear regressor for the emulation tasks of the 2nd 
order system (upper diagram), 10th order system (middle diagram), and 
Volterra task (lower diagram). In the plots, for the 2nd and 10th order 
systems, and the Voletrra task, the MSEs of the system output were 
1.89 X 10-', 3.80 X 10-^, and 8.90 x lO"'', respectively and those of the 
linear regressor were 9.96 xlO"', 4.95 xlO-**, and 9.48 xlO""*, 
respectively. For each task, the system output showed better performance 
than the linear regressor. Note that the output of the linear regressor is not 
a straight line, but a scaled version of the input with an offset (this outcome 
is due to the scaling of the figure). 



MSEs over 20 trials. In each task, our system showed significantly 
low MSE. These results suggest that our system is able to exploit 
the non-linearity and memory originates fi-om the passive body 
dynamics of the soft robotic arm to perform the task. 

Unlike a conventional computational network or device, our 
system receives input as a mechanical rotation of the base, which 
generates the physical motion of the arm. Thus, we can easily 
imagine that, for example, if the input scaling of the base rotation 
range, R, is small, then the arm can only vibrate slightly, which 
does not generate diverse body dynamics. Moreover, since the 



Table 2 | Comparisons of MSE between the output of the system and 
that of the linear regressor for each task. 





System output 


Linear regressor 


p value 


2nd order 


1.84 ±0.05 


1.01 ± 0.03 


p < 0.001 




(xlO-') 






10th order 


3.77 ± 0.03 


5.05 ±0.11 


p < 0.001 




(x10-6) 


(xlO-'*) 




Volterra task 


8.89 ±0.18 


9.28 ±0.11 


p < 0.001 




(x10-i^) 


(xlO-4) 





input is linearly scaled to the base rotation range in our setting, 
the degree of this scaling changes not only the amplitude of the 
rotation but also the speed of the rotation. Considering that the 
water friction on the arm shows non-linear dependence on the 
velocity and the angle of the compartments (Equations 12, 13, 
and 14), the property of the body dynamics would change accord- 
ing to the degree of input scaling, R. Accordingly, the performance 
of our system would also change for each task. In order to validate 
this, we varied R from 15 to 90, and observed how the perfor- 
mance of the system changed for each R. Figure 5A shows the 
results. We can confirm the different individual error profile with 
respect to R for each task. First, small R values (around 15) show 
the highest errors; errors gradually start to decrease according to 
increases in R values in all tasks. But, for example, in the case of 
the 2nd order system, the error starts to increase again at around 
R = 30 and has a local maximum at around R = 45. In the case 
of the 10th order system, the error just decreases monotonically 
as the value of R increases. In the case of the Volterra task, the 
error shows the minimum at around R= 55 and start to increase 
monotonically as the value of R increases. This suggests that, even 
if the mechanical structure of the arm is the same, certain behav- 
iors of the arm can reveal especially high computational power in 
some tasks, but not in others. 

As explained in section 1, our system is essentially classi- 
fied as a reservoir computing approach, in which a number of 
randomly coupled nodes are usually used as a computational 
resource, and where each node has a statistically uniform role. 
On the other hand, in our system, due to the intrinsic body struc- 
ture, we can expect that there is a specific role for each body 
part. Here, we aim to investigate this point in two ways. First, 
when running the evaluation phase, we take out the readouts 
from one compartment and analyze the error. Note that, in this 
analysis, we use the readouts, which are trained with 20 com- 
partments. By iterating this procedure for each compartment, 
we can investigate how each compartment contributed to the 
task performance when the readouts were fully connected. We 
call this contribution ratio analysis of the compartments. Second, 
we perform the entire experiment (i.e., washout, learning and 
evaluation phase) with only 19 compartments by skipping one 
compartment, and compare the performance with that obtained 
using 20 compartments. The difference with the previous contri- 
bution ratio analysis is that the readouts are, in the first place, 
trained to maximize the performance with 19 compartments 
including the entirely new readout weights. Since the readout 
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FIGURE 5 I (A) Dependence of the performance of the system output (MSE) 
on the base rotation range, R, for each task. The plots show the averaged 
values of MSE over 20 trials. The error bars show the standard deviation. The 
averaged value of MSE of the linear regressor was much higher than that of 
the system output with respect to each R value (usually an order of 
magnitude higher). (B) Results of the contribution ratio analysis for each task. 
The horizontal axis shows the number of the compartment discarded in the 



evaluation phase, and the vertical axis shows the contribution ratio (CR). The 
plots show the averaged CR values over 20 trials, and the error bars show the 
standard deviations. (C) Results of the computational power analysis for each 
task. The horizontal axis shows the number of the compartment excluded to 
train the readout weights, and the vertical axis shows the performance ratio 
(PR). The plots show the averaged PR values over 20 trials, and the error bars 
show the standard deviations. See text for the details. 



weights are optimized without using a specific compartment, we 
can infer back the overaO computational power of each com- 
partment as a deviation from the original 20 compartment case. 
We call this computational power analysis of the compartments. 
Hereafter, the base rotation range, R, is fixed to 60 for the 
analyses. 

In the contribution ratio analysis, the experimental procedure 
is the same as that explained in section 2.2, except that the eval- 
uation phase is performed by taking out the readouts from a 
specific compartment (thus, four nodes are excluded). We iter- 
ate this procedure for each compartment by using the same input 
time series and body dynamics in a trial and calculate the MSE 
for each case. After testing all the compartments, we normalize 
them with the maximum MSE collected, and obtain contribu- 
tion ratio (CR) for each compartment. If the CR is high for the 
compartment, then it implies that this compartment was con- 
tributing to the task performance largely when the readouts were 
fully connected. Figure 5B shows the result of the contribution 
ratio analysis for each task over 20 trials. In the case of the 2nd 



order system, although compartments 2, 3, 4, 8, and 9 seem to 
have high CRs, the standard deviations for these are also high, 
while compartments 1, 6, 19, and 20 have low CRs with low stan- 
dard deviations. This suggests that specific compartments, such 
as 1, 6, 19, and 20 always contribute less to the task performance, 
while the computational role for this task is relatively distributed 
throughout the resting compartments, and among them, there is 
no specific compartment that consistently has a high contribu- 
tion. In the case of the 10th order system and the Volterra task, 
the situation is different. There seems to be key compartments 
that always show high contributions to the task performance. 
In the case of the 10th order system, compartments 16, 17, and 
18 show high performance, and in the case of the Volterra sys- 
tem, compartments 14 and 15 show high performance. Overall, 
these results suggest that our system adopts various strategies in 
the performance of computational abilities, according to the task. 
One strategy is to distribute the computational role throughout 
the entire arm, while the other is to always select and rely on the 
motion of specific body parts. 
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Next, in the computational power analysis, the experiment is 
performed both under the default condition (using 20 compart- 
ments) and with the exclusion of a specific compartment when 
training the readout weights. This is done by using the same 
input time series and body dynamics in a trial for each task. We 
calculated the MSEs in the evaluation phase for both cases and 
divided the MSE of the case without the specific compartment 
by that of the default condition and obtained the performance 
ratio (PR) for each task in each trial. Thus, if the PR is larger 
than 1.0, it implies that the task performance is worse than in 
the default condition, and the value indicates the degree of how 
the exclusion of the specific compartment affects the overall task 
performance in terms of ratio. Figure 5C shows the result of the 
computational power analysis for each task over 20 trials. In the 
2nd order system task, we can clearly see that there are high 
PR values around the base of the arm (in compartments 1, 2, 
3, and 4) suggesting that these compartments contain signifi- 
cant information for the performance of this task. Similarly, in 
the Volterra task, there are high PR values in compartments 14, 
15, and 16. On the other hand, in the 10th order system task, 
PR values lower than 1.0 are shown around the tip of the arm 
(in compartment numbers higher than 15), suggesting that these 
compartments have a negative influence on the performance of 
this task. We speculate that this is caused by an overfitting effect 
produced by the compartments around the tip of the arm. The 
network is too specialized for the learning data and is not able 
to generalize to the new (evaluation) data. The reason for this 
should be explored in more detail in future work. Overall, we 
showed that there are specific regions in the body parts that con- 
tain positive or negative information for the performance of the 
tasks. 

In this section, we first demonstrated that our soft robotic 
arm can perform the task to emulate non-linear dynamic systems 
by positively exploiting the non-linearity and memory originates 
from its body dynamics. We also confirmed that the way the input 
is applied (in our case, the amplitude range for the arm move- 
ment) significantly affects the computational ability, and its body 
parts show specialized roles due to their intrinsic morphologi- 
cal structure and corresponding diverse body dynamics, unlike 
the conventional reservoir computing approach. In the next sec- 
tion, we see how these body dynamics can potentially be used to 
control the arm's motion in a closed-loop manner by embedding 
non-linear limit cycles. 

3.2. TASK 2: CLOSED-LOOP CONTROL— EMBEDDING NON-LINEAR 
LIMIT CYCLES 

In this section, we show the results for Task 2. By following the 
procedure described in section 2.2.2, we conducted a number 
of computer simulations to train the readouts with various val- 
ues of the base rotation range R and the degree of white noise 
V. As a result, we heuristicaUy found that the system perfor- 
mance is extremely sensitive to the setting of these parameters [as 
opposed to the results presented for the simpler and abstract net- 
works used in Hauser et al. (2012)]. (As for R, we have already 
shown in the previous section that R changes the computational 
power of the system significantly.) If these parameters were not 
set appropriately, we often observed that, when the system was 



switched from the teacher forcing condition to the closed-loop 
control, the arm gradually approached the resting state or showed 
unrealistic behaviors due to numerical problems. For the latter 
case, since we adopt the position control of the base angle, if 
the output showed much higher values than that of one timestep 
before (for example, if |0(f-|- 1) - 0(0 1 > 10, then the arm 
would have to rotate its base extremely quickly, namely, larger 
than lO^deg/s, which is unrealistic in the physical platform), 
then, as a result, the simulator showed numerical problems. We 
carefully discarded these cases from our experiment. Even if the 
system has a high computational power as we saw in the previous 
section, the closed-loop setting requires additional care due to the 
stability issues. Since the output, which includes the error, is fed 
back to the system as input, the error may grow larger and larger 
in each simulation timestep. 

Figures 6-8 show the typical results we obtained when the 
arm does not approach the resting state or the unrealistic behav- 
iors mentioned above, for closed-loop control of the Van der Pol 
limit cycle, the quadratic limit cycle, and the Lissajous curve. 




system output 




output 



Xi(t)/0,(t) 

FIGURE 6 I Results for implementing the Van der Pol limit cycle. At 

timestep 5000, tlie system is switclied from the teacher forcing condition 
to the closed-loop control (black line). (A) The time series of the lengths of 
springs 1, 2, 3, and 4 for all 20 compartments. (B) Comparisons between 
the system output (Oi (t) and 02(t) (blue lines)) and the target output (xi (t) 
and X2(t) (red lines)). (C) Comparisons between the system output (blue 
lines) and the target output (red lines) in the Oi (f) - 02(f) (xi (f) - X2(t)) 
plane. The time series from timestep 5000 to timestep 70,000 are overiaid. 
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respectively. '3 For the parameters (_R, v), we adopted (_R, v) = 
(130, 1.5 X 10"*^), (90, 1.0 X 10"*'), and (10, 1.0 x 10""), for 
the Van der Pol limit cycle, the quadratic limit cycle, and the 
Lissajous curve, respectively. For the Van der Pol limit cycle, 
we can see that the system is not implementing the target tra- 
jectory (Figures 6B,C), but is rather implementing an irregular 
one (Figure 6C). We also observed that the behavior of this 
trajectory is not as stable, but rather constantly changes its 
trajectory for each cycle, and this change remains throughout 
the trial. However, the results for the quadratic limit cycle and 
the Lissajous curve show almost a complete fit with the target 
trajectory (Figures 7B,C, and 8B,C, respectively). We confirmed 



^Strictly speaking, to prove that the embedded trajectory is really a "limit 
cycle", we need to analytically show whether the trajectory is an attractor of 
the system. In our case, this is unrealistic because the equation governing the 
mechanical system is too complex, and we would have to rely on a heuristic 
approach. Accordingly, as we see later, we here call the embedded trajectory a 
"limit cycle" if and only if the trajectory can stay at the target limit cycle for 
1,000,000 timesteps and the trajectory has a certain attraction when perturbed 
externally. 




C system output 




output 



-TOO 0 100 



Xt(t)/0,(t) 

FIGURE 7 I Results for implementing the quadratic limit cycle. At 

timestep 5000, the system is switclied from the teacher forcing condition 
to the closed-loop control (black line). (A) The time series of the lengths of 
springs 1,2,3, and 4 for all 20 compartments. (B) Comparisons between 
the system output (Oi (t) and 02(t) (blue lines)) and the target output (xi (t) 
and X2(t) (red lines)). (C) Comparisons between the system output (blue 
lines) and the target output (red lines) in the Oi (f) - 02(f) (xi (t) - X2(t)) 
plane. The time series from timestep 5000 to timestep 70,000 are overlaid. 



that these trajectories were stable enough to run for 1,000,000 
timesteps without leaving the trajectories of the target limit cycles. 
These results suggest that the task performance of the closed-loop 
control is not only restricted to the degree of non-linearity or 
memory required for the limit cycles but is also dependent on how 
the arm is driven. These preferences are caused by the intrinsic 
structure of the body. From now on, by using the system embed- 
ding the quadratic limit cycle (Figure 7) and the Lissajous curve 
(Figure 8), we move on to analyze the stability of the closed-loop 
controls and the role of each body part in these limit cycles as we 
saw in the previous section. 

One important aspect to evaluate the embedded closed-loop 
control is its robustness against external perturbations. To test 
this, we added white noise, 80i(f) and 802(f)) in the range of 
€ (80i,2(f) e[— €, €]) to the two motor outputs of the embed- 
ded control, such as Oi(t) -|- 8Oi(0 and 02(f) + 802(f), during 
t = 6000-7000 timesteps for each as an example. Figure 9 shows 
the typical results of the performance of the embedded quadratic 
limit cycle regarding each noise level (e = 10.0, 5.0, and 1.0). 
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FIGURE 8 I Results for implementing the Lissajous curve. At timestep 
5000, the system is switched from the teacher forcing condition to the 
closed-loop control (black line). (A). The time series of the lengths of 
springs 1, 2, 3, and 4 for all 20 compartments. (B) Comparisons between 
the system output (Oi (t) and 02(0 (blue lines)) and the target output 
(xt (0 and X2(0 (red lines)). (C) Comparisons between the system output 
(blue lines) and the target output (red lines) in the Oi (0 - 02(f) (xi (f) - X2(f)) 
plane. The time series from timestep 5000 to timestep 70,000 are overiaid. 
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FIGURE 9 I Typical examples showing the performance of the embedded 
quadratic limit cycle with external noise. The white noises, SOi (t) and 
502(t), in the range of c (50i_2(f) el— e, el) are added to the motor outputs 
during f = 6000 - 7000 timesteps. (A) Typical trajectories of Oi (t) and 02(0 



compared with the target time series. The upper, middle, and lower 
figures show the case with e = 10.0, 5.0, and 1.0, respectively. (B) The 
corresponding plots in the Oi (f) - 02(f) (xi (t) - X2(t)) plane. For each noise 
level, the trajectories are overlaid from 5000 to 30,000 timesteps. 



We can clearly see that even if the noise level is relatively large, 
such as € = 10.0, the trajectories eventually recover to the limit 
cycle, which suggests that the embedded quadratic limit cycle is 
robust against external noise. We also confirmed that, even if we 
elongate the duration of time for the added noise, the system can 
successfully recover its performance. Note that, although the per- 
turbed trajectories came back toward the limit cycle (Figure 9B), 
the oscillation phase was often shifted (Figure 9A). This was 
mainly caused by the relatively long duration of time for adding 
noise. We observed that by shortening this duration, this phase 
shift tendency can be reduced accordingly. Next, let us see the case 
for the embedded Lissajous curve. Compared to the quadratic 
limit cycle case, the system is less robust. When € is more than 
around 0.3, we often observed that the trajectories go out from 
the limit cycle and never come back. Figure 10 shows the typi- 
cal results of the performance of the embedded Lissajous curve 
for each noise level less than 0.3 (e = 0.3, 0.15, and 0.1). If the 
noise level was less than around 0.3, we observed the system per- 
formance recovered toward the limit cycle as in the quadratic 



limit cycle case. However, even in this noise range, we some- 
times observed an unstable trajectory as shown in Figure 11. In 
addition, similarly to the quadratic limit cycle case, even if the per- 
turbed trajectories came back toward the limit cycle (Figure lOB), 
the oscillation phase was often shifted (Figure lOA). 

Now, we move on to see the role of each body part (compart- 
ments or springs) as we saw in Task 1. Since the motor outputs 
and the body dynamics are reciprocally coupled through the feed- 
back loop, the scheme we adopted in the Task 1 case, such as 
skipping one compartment, will cause unrealistic behavior of the 
arm due to numerical problems, as explained previously, and can- 
not always be adopted to appropriately evaluate the contributions 
of the body parts. Accordingly, we aim to investigate the contri- 
bution of each body part in terms of robustness against noise. 
Namely, we evaluate how the noise added to each body part affects 
the overall system performance. As is obvious from the system 
construction (Equations 24, 25, and 26) for the closed-loop con- 
trol, the slight difference in the motor outputs at timestep f can 
affect the corresponding sensory time series, i.e., the lengths of the 
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outputs during f = 6000 — 7000 timesteps. (A) Typical trajectories of Oi (f) level, the trajectories are overlaid from 5000 to 30,000 timesteps. 



springs, and this effect influences the motor outputs at timestep 
f + 1. That is, according to the slight difference in the motor out- 
puts, 80i(f) and802(f)> the sensory time series, Sy(f), will deviate 
from the original expressed as s'--(t) = 5y(f) -|- 8sy(f), where s'-j(t) 
is the actual spring length at timestep t. Then, from Equation (26) , 
the outputs at timestep t + I can be simply expressed as: 

20 4 

0'(f+l) = I]I]i^ou.4(f), (30) 

;=1,=1 
20 4 

= I]I]^ou,(^.;(f) + S5y(f)), (31) 

i=l,=l 

20 4 

= 0(t+l) + J2Yl i^outSsy(f), (32) 

] = U=\ 

= 0(f-M)-F80(f-|- 1), (33) 



where 0'(f) and 0(f) are the actual and original motor outputs at 
timestep f, respectively. Note that for simplicity we dropped the 
index expressing two outputs. Since the deviation of the motor 

outputs is expressed as 80(t-|- 1) = J2j=i X!f=i ^out^^yCO. we 
can investigate how the noise applied to a single spring at timestep 
t can affect the motor outputs at timestep f -|- 1 by fixing the other 
springs as the original. By investigating how 80(f -|- 1) evolves 
through time, we can also evaluate the effect of the noise against 
the overall system performance. 

Now, let us assume that the noise was applied to the sen- 
sory value of the spring i in compartment ; at timestep f by 
fixing the other sensory values as the original. Then, the deviation 
of the motor output at timestep t + I can be simply expressed 
as 80(f -|- 1) = w'^^^?>Sij{t), which straight-forwardly means that 
the degree of 80(f -|- 1) is only linearly dependent on the read- 
out weight of the focused spring. Therefore, we can infer and 
compare how the noise added to each sensory value affects the 
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FIGURE 11 I Unstable trajectory of the embedded Lissajous curve with 
external noise in the range of 6 = 0.2. The plot is in tlie Oi (f) - 02(f) 
(xi (f) - X2(t)) plane. The trajectories are overlaid from timestep 5000 to 
tinnestep 30,000. The target trajectory is also shown as a reference. 



motor output at timestep f + 1, regarding the fixed noise value, 
only by checking the weight distributions. Simply saying, if the 
value of IWoutI is large, then the effect of the noise for spring 
i in compartment j at timestep t on the outputs for timestep 
t + 1 is also large, which means this spring makes a big con- 
tribution to the transition of the outputs for timestep f + 1. 
Figure 12 shows the readout weight distributions of Wout.i and 
Wout,2 for the embedded quadratic limit cycle (Figure 12A) and 
Lissajous curve (Figure 12B). We can see that the distributions 
show a characteristic pattern for each limit cycle. For exam- 
ple, in the embedded quadratic limit cycle case, the value of 
the weights often seems to have a symmetric and correspond- 
ing distribution over springs in each compartment (e.g., between 
spring 1,2 and spring 3,4 in compartment 3-13), and even over 
Wout.i and Wout,2 (Figure 12A). In the embedded Lissajous curve 
case, this type of symmetric and corresponding distribution 
can also be found within each weight, but not over Wout.i and 
Wout,2 (Figure 12B). In addition, for the distribution of Wout.i, 
the value of the weights is almost zero in compartment 7-20, 
which means the external noise applied to these compartments 
at timestep t will not affect the motor outputs at timestep t + I 
very much. 

To systematically proceed with this line of analysis for the 
overall system performance, we need to confirm whether the 
large deviation in the sensory value at a specific time also leads 
to the large deviation of motor outputs over time. Although this 
seems to be trivial in our system, it is not trivial in general because 
the transition 80(f) 8sy(f) depends on the construction of 
the body (for example, imagine the sensory value that exhibits 



a saturation)''. To evaluate this, we added a small white noise in 
the degree of € to the motor outputs only at timestep 6000, and 
investigated how the differences in the motor outputs, |80(f)l, 
and the spring lengths, |8sy(t)|, carry on over time according to 
the degree of € by measuring the mean square errors between 
the actual and original motor commands, such as MSEqi = 

T Ef=i(o'i(f) - Oi(f))^ MSE02 = ^ ELi(02(f) - 02(f))', 

and between the actual and original sensory time series as 
MSEspring = Ell E ■=! Ell (4(0 - s^j{t))\ where 
T = 500 throughout this experiment. Note that we consider only 
a small range of noise around 6 e [0.005, 0.1], since if the noise 
level is too large, the trajectories often show phase shifts as we saw 
in the previous analysis (Figures 9, 10), which make the measures 
miss capturing the intended difference even if the trajectories 
were in the original limit cycles. Figure 13A shows the results of 
the averaged MSEqi, MSE02, and MSEsp^ing for each embedded 
limit cycle. We can clearly confirm that according to the increase 
in the noise level €, the value of each measure also increases, 
which means that the large deviation of the motor outputs at a 
specific time also leads to a large deviation in the motor outputs 
over time. In addition, in the embedded Lissajous curve case, 
MSEoi is larger than MSE02 for each t value, which suggests 
that the output Oi(f) is more sensitive than 02(f) (Figure 13A 
(right)). 

Now, we are ready to investigate the role of each body 
part. According to the results shown in Figure 13A, the weights 
assigned for each spring directly reflect how the noise added 
to each sensory value affects the overall system performance. 
To correspond to the results in Figure 13A, we first calculated 

dw = ^(M'out 1)' + (''^out 2)' ^'"-h spring i in each compart- 
ment j, since the size of the noise to the motor output can be 

expressed as 8syy(w;^^t i)2 + (wl^^_^)^ = ySOKfPTSO^ < 
€ by scaling the size of 8sy in the appropriate range. Thus, the 
value of dw directly reflects the contribution of each spring to 
the overall system performance regarding the fixed noise value. 
Figure 13B shows the value of dw for each spring i accord- 
ing to each compartment for the embedded quadratic limit 
cycle and Lissajous curve. Interestingly, the value of dw for 
each spring is almost the same within each compartment for 
both limit cycles, which means that the contribution of each 
body part can be expressed at the compartment level. In the 
case of the quadratic limit cycle, the value of dw shows almost 
a zig-zag pattern and gradually decreases when the compart- 
ment number increases from the base toward the tip [Figure 13B 
(left)]. In the case of the Lissajous curve, only the compart- 
ment around the base shows high values for dw [Figure 13B 
(right)]. These results suggest that, according to the limit cycles 
embedded, the sensitivity and degree to affect the overall system 
performance against the noise show different tendencies for each 
body part. 



^Theoretically, this is to analyze the basin structure surrounding the original 
trajectory. Due to the number of parameters, instead of analyzing the basin 
structure according to 80i(t) and 502(t), we analyzed the basin volume in 
terms of the error measures introduced later according to |SO(t)|. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



July 2013 I Volume 7 1 Article 91 | 15 



Nakajima et al. 



A soft body as a reservoir 



A W 



out. 



spring 

12 3 4 














































1 




































1 









































































1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 

^ , „ compartment 



Q-C\l 



I 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 

compartment 

B 






1 


2 

lilt 


3 4 5 6 
? 


7 


8 9 10 11 12 13 14 15 16 17 18 1920 
compartment 


■<ij- 
























I 




























































i 

CLCM 
W 




























■ 


















■ 








































1 


2 


3 


4 


5 


6 


7 


8 9 10 11 12 13 14 15 16 17 18 1920 
compartment 



XI 04 

4.0 
2.0 
0 

-2.0 
1-4.0 

xlO" 
I 4.0 

2.0 

0 

-2.0 
-4.0 

X108 
2.0 
1.0 
0 

-1.0 
-2.0 
-3.0 

X 106 

1.5 



- 0 



"-1.5 



FIGURE 12 I (A) Plots showing tine readout weiglit distribution for the quadratic limit cycle. (B) Plots showing the readout weight distribution for the Lissajous 
curve. For each figure, the upper and lower graphs show w„u, i and w„„, 2, respectively. 



In this section, we have investigated whether our soft robotic 
arm can embed Kmit cycles in a closed-loop manner, and 
have shown that several properties, such as the system perfor- 
mance, the robustness against external noise, and the role of 
each body part differ according to which limit cycle to embed. 
Since we are adjusting only the linear readouts by using the 
same body, we can speculate that these specific properties cor- 
responding to each limit cycle are caused by the intrinsic body 
structure. 

4. DISCUSSION 

In this paper, by using the dynamic simulator of the soft robotic 
arm inspired by the octopus, we demonstrated that the robot's 
body dynamics are already capable of emulating non-linear 
dynamical systems and embedding non-linear limit cycles in 
a closed-loop manner by only adjusting the fixed linear read- 
outs. The arm we used did not contain any rigid components. 
Instead, it is soft, including only springs, which are aligned to 
mimic the muscular structure of the octopus. This resulted in 



several compartments, each of which had a specific muscular- 
hydrostat property, which enforced the springs to be coupled 
in well-defined, but constrained, manner. In addition, the arm 
was assumed to be immersed in an underwater environment, 
in which the friction constants were identified via CFD simu- 
lations. All these factors, including this intrinsic body structure 
and its interaction with the environment, generated diverse body 
dynamics, including rich non-linearity and memory. The tech- 
nique presented here allowed us to exploit these properties as 
computational resources. In addition, it is possible to infer the 
amount of non-linearity and memory that can be potentially 
exploited for information processing in terms of the task perfor- 
mance. For roboticists, this may open up the way to quantitatively 
characterize which control is efficient for which body design, as 
well as outsourcing the control load to the body parts. Although 
we kept the arm's mechanical structure as bio-inspired as possible 
throughout the analyses, it would also be meaningful to inves- 
tigate how the information processing capability would change 
if the arm's mechanical properties (such as stiffness, damping. 
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drag parameters) are altered. This line of experimentation will be 
included in our future work. 

For the emulation tasks of non-linear dynamical systems, 
in addition to its high computational power, we showed that 
each body part has a specific role according to the task type. 
Additionally, for the closed-loop control tasks, we showed that 
the arm prefers some limit cycles over others (i.e., the quadratic 
limit cycle and the LIssajous curve were possible to embed, while 
the Van der Pol limit cycle was not). These obvious and specific 
coherencies are usually not observed in conventional reservoir 
computing where the reservoir consists of randomly coupled 
non-linear computational elements, suggesting that these prop- 
erties originate from the intrinsic body structure. 



From a biological systems point of view, this result seems 
natural. In nature, animals adapt to their respective ecologi- 
cal niches, where they evolve their body morphology to survive 
within their environment. The octopus is not an exception; its 
specific body structure is specialized to permit survival in a com- 
plicated underwater environment, enabling it to behave efficiently 
in particular ways. In this context, it would be interesting to 
investigate whether the arm could embed more biologically plau- 
sible behaviors in future work. For example, as we mentioned 
earlier, it is well known that the octopus adopts a specific strat- 
egy for reaching, called bend propagation (Gutfreund et al., 1996; 
Gutfreund, 1998; Sumbre et al., 2001; Yekutieli et al., 2005a,b). 
In this specific motion, it is suggested that the CNS only ini- 
tiates the motion and all the muscle activations are handled at 
the PNS level (Gutfreund et al, 1996; Gutfreund, 1998; Sumbre 
et al., 2001). Several researches have investigated this behavior by 
directly extracting the muscle contraction patterns from the real 
octopus, and by externally applying these patterns to the octopus 
arm models (Gutfreund et al., 1996; Gutfreund, 1998; Yekutieli 
et al., 2005a,b). On this point, our technique presented here may 
reveal further insights on this overall scheme by including the role 
of the arm's body dynamics. Considering that the PNS does not 
have a plasticity (Kandel et al, 2000), it would be worth inves- 
tigating how the arm's body dynamics, together with the PNS, 
modeled as a linear and static feedback loop onto the arm, embeds 
the motor patterns of bend propagation according to the initia- 
tion command sent by the CNS. This line of experiment can be 
investigated in future work. 

There exists a growing number of documented cases in nature, 
which support that certain morphologies found in animals are 
facilitating a kind of computation. This observation is usually 
characterized by the term morphological computation. For exam- 
ple, the non-linear, non-homogeneous spatial arrangement of the 
ommatidia in insect eyes are more dense toward the front than 
on the side in order to compensate for motion parallax, which 
is non-linear (Franceschini et al, 1992). The morphology coun- 
teracts the non-linearity introduced by the parallax; hence, the 
complexity of the computational tasks to steer through obstacles 
based on the visual input is reduced. Since the resulting task for 
the brain is now simpler and not non-linear anymore due to the 
"clever" morphology, one could argue that part of the computa- 
tion is conducted by the morphology. While this is a very simple 
case of a morphological computation, since the given morphol- 
ogy represents only a static, non-linear mapping, the concept does 
go further, if we consider, for example, soft, compliant bodies. 
Such bodies exhibit interesting dynamic properties, such as fading 
memory and non-linearity. Examples of such complex computa- 
tions outsourced to the physical layer are passive walkers (Collins 
et al., 2005). Their design pushes the limits of what can be out- 
sourced to the physical body, in so far that no controller (i.e., 
CPU) is needed at all. The mechanical design inspired by the mus- 
culoskeletal structure enabling "preflexes", which can self-stabilize 
movements through its elastic material properties, also gives such 
an example (Brown et al., 1995; Blickhan et al, 2007; Proctor 
and Holmes, 2010). Their morphology (i.e., the mechanical, soft 
design and the environment) is able to "do" all the computations 
needed to walk robustly. While such robots are impressive, their 
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disadvantage is their inflexibility by restricting the computation 
to a fixed physical body only. In a biological system, a sensible 
distribution of the computation between the body and the brain 
is more probable. Despite the number of biological examples and 
the series of robots, which have been built considering the concept 
of morphological computation in their design, there are still a few 
studies characterizing the concept within a quantitative frame- 
work (Hauser et al, 2011, 2012; Fuchslin et al., 2013). In this 
context, we believe that the approach presented here would be 
one of the interesting directions for further study. 
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